ELSEVIER 


Available  online  at  www.sciencedirect.com 


SCIENCE 


DIRECT® 


Journal  of  Power  Sources  157  (2006)  630-640 


www.elsevier.com /locate /jpowsour 


Multi-level  reduced-order  thermal  modeling  of  electrochemical  capacitors 

Ph.  Guillemet*,  Y.  Scudeller,  Th.  Brousse 

Laboratoire  de  Genie  des  Materiaux,  Ecole  Polytechnique  de  I’Universite  de  Nantes,  La  Chantrerie,  rue  Christian  Pauc,  BP  50609,  44306  Nantes  Cedex  3,  France 

Received  7  February  2005;  received  in  revised  form  14  July  2005;  accepted  20  July  2005 

Available  online  17  October  2005 


Abstract 

Ultra  capacitors  are  major  components  used  as  power  sources  offering  a  combination  of  high  power  and  high  energy.  Thermal  management  of 
ultra  capacitors  is  one  of  the  main  issues  for  the  design  of  safe  powerful  systems.  This  paper  presents  multi-level  electrothermal  modeling  that 
can  be  used  to  design  ultra  capacitor  structures  meeting  reliability  requirements  of  power  applications.  The  multi-level  modeling  is  based  on  both 
numerical  and  analytical  approaches  enabling  us  to  take  into  account  different  scales  through  finite  element  method  computations,  shell  network 
models,  homogenization  methods  and  ultra-reduced-order  model.  Basic  understanding  of  electrothermal  behavior  is  performed.  Influence  of  cell 
characteristics  and  cooling  conditions  are  studied. 

©  2005  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Electrochemical  capacitors,  so-called  supercapacitors  or  ultra 
capacitors,  are  reversible  energy  storage  devices  intermedi¬ 
ate  between  secondary  batteries  and  conventional  capacitors 
(Conway  [1]).  Basically,  they  offer  instantaneous  power  den¬ 
sity  more  significant  than  that  of  batteries  and  energy  density 
larger  than  that  of  dielectric  capacitors.  Additionally,  they  can 
be  charged/discharged  more  than  105  times  without  significant 
energy  loss.  Ultra  capacitors  are  attractive  for  a  wide  variety 
of  power  applications  such  as  telecommunication  satellites  or 
hybrid  electrical  vehicles.  In  these  systems,  pulsed  power  con¬ 
version  and  long-term  charge-discharge  cycling  are  required  in 
combination  with  standard  secondary  batteries. 

Three  main  classes  of  materials  are  used  to  prepare  the 
electrodes  of  electrochemical  capacitors:  carbon  (Becker  [2]), 
electronically  conducting  polymers  (Gottesfeld  et  al.  [3])  and 
metal  oxides  such  as  RuC>2  (Trasatti  and  Buzzanca  [4])  or 
more  recently  MnC>2  (Lee  and  Goodenough  [5]).  These  elec¬ 
trode  materials  are  used  in  aqueous  based  electrolytes  (H2SO4, 
K2SO4,  etc.)  or  non-aqueous  electrolytes  that  are  character¬ 
ized  by  a  wider  electrochemical  stability  window  relative  to 
aqueous  electrolytes.  Thus,  the  combination  of  a  solvent  such 
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as  acetonitrile  or  propylene  carbonate  with  different  salts  (e.g. 
tetraethylammonium  tetrafluoroborate)  enabled  the  increase  of 
the  maximum  cell  potential  up  to  3  V  (Conway  [1]).  However, 
despite  the  remarkable  performance  of  these  organic-based  sys¬ 
tems,  they  suffer  from  the  use  of  highly  toxic  and/or  flammable 
solvents  which  can  cause  severe  safety  hazards.  Subsequently, 
thermal  management  of  the  supercapacitors  is  the  main  issue  for 
the  design  of  safe  powerful  systems.  This  management  is  actu¬ 
ally  barely  conducted  by  commercial  software  that  do  not  take 
into  account  the  internal  layer  structure  of  ultra  capacitors. 

Pulse  power  has  durations  ranging  from  about  10-3  s  to  100  s. 
The  move  to  higher  power  will  continue  in  the  future,  and  it  is 
now  desirable  to  develop  a  fundamental  understanding  of  elec¬ 
trothermal  behavior  of  ultra  capacitors. 

The  general  aim  of  this  study  is  to  develop  multi-level 
reduced-order  thermal  modeling  of  ultra  capacitors,  including 
both  material  properties,  structure  of  electrochemical  cells  and 
packaging.  The  objective,  besides  a  fundamental  description  of 
ultra  capacitor  thermal  behavior,  is  to  bring  up  a  useful  tool  for 
the  design  of  ultra  capacitors  and  to  assess  cooling  conditions. 

2.  Temperature  and  heat  production  effects  on 
performance 

Power  losses  are  produced  by  charge-discharge  current 
cycles  which  cause  undesirable  hot  spots  affecting  both  reliabil- 
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Nomenclature 

c  specific  heat  (J  K- 1  kg- 1 ) 

d  distance  (m) 

e  thickness  (m) 

h  heat  transfer  coefficient  (W  m-2  K- 1 ) 

Kj  terminal  thermal  conductance  (W  K-1  m-2) 
/  distance  (m) 

N  number  of  cells 

Qtot  total  heat  source  (W) 

Qy  volume  heat  generation  rate  (W  m-3) 

R  thermal  resistance  (K  W- 1 ) 

t  time 

T  temperature  (K) 

v,  y,  z  space  variables 

Greek  letters 

a ,  dimensionless  coefficients 

A  thermal  conductivity  (W  m-1  K_1) 

p  density  (kg  m_ 3 ) 

0  heat  flux  (W) 

Subscripts 
amb  ambient 

c  collector 

e  electrode 

e+  positive  electrode 

e—  negative  electrode 

es  electrode- separator 

i  index 

is  insulator 

s  surface  package 

tot  total 

w  terminal  wire 


ity  and  performance.  Besides  ageing  of  supercapacitor  elements, 
overheating  also  induces  electrochemical  variations  of  internal 
characteristics. 

Temperature  influences  greatly  the  electrochemical  behav¬ 
ior  of  components,  like  self-discharge.  The  leakage  resistance 
being  higher  than  the  one  of  conventional  capacitors,  thermal 
degradation  may  occur.  Capacitance  and  leakage  resistance  of 
individual  cells  have  to  be  very  constant  during  the  whole  life 
of  the  device.  A  uniform  temperature  is  therefore  desirable. 

Moreover,  ultra  capacitors  have  the  disadvantage  of  rather 
low  voltage  (1-3  V),  so  many  single  cells  have  to  be  connected  in 
series  in  order  to  achieve  several  hundred  volt  capacitor  voltage. 
Capacitors  can  be  destroyed  if  exposed  to  transient  over  voltages. 
In  a  series  connected  string,  some  cells  may  be  submitted  to  an 
over  voltage  and  thermal  damage. 

Electrothermal  behavior  of  ultra  capacitors  can  be  diffi¬ 
cult  to  predict  because  a  large  series  of  transport  phenomena 
(ionic  and  electronic  transport,  heat  and  mass  diffusion)  and 
structure  heterogeneities  are  involved.  Ultra  capacitors  form 
tri-dimensional  non-isotropic  structures  basically  composed  of 
interconnected  arrays  of  electrochemical  unit  cells  stacking  up 


Electrodes  :  actived 


to  several  hundred  cylindrical  or  prismatic  unit  cells.  A  wide 
range  of  materials  is  used  (composite  electrodes,  electrolyte  and 
porous  separator,  metal  current  collectors,  interconnects,  pack¬ 
age).  Temperature  profiles  depend  upon  the  internally  generated 
power  losses  that  can  be  non-uniform  within  electrochemical 
cells.  Current  profiles  frequency,  magnitude  of  current,  unit  cell 
structure,  material  properties,  arrangement  of  cells,  packaging 
and  cooling  conditions  have  a  great  influence  on  the  working 
temperature. 

3.  Description  of  the  electrochemical  unit  cells 

Ultra  capacitor  bipolar  cells  consist  of  high  surface  area  elec¬ 
trode  materials  loaded  with  electrolyte  (Fig.  1).  Negative  and 
positive  electrodes  are  separated  by  a  porous  separator  impreg¬ 
nated  with  liquid  electrolyte. 

Thick  film  electrodes  are  composite  materials  typically  pre¬ 
pared  by  mixing  80%  active  powder,  7.5%  graphite,  7.5%  acety¬ 
lene  black  and  5%  polymer  (PTFE). 

Symmetric  capacitors  are  assembled  from  activated  carbon 
powder.  Hybrid  capacitors  are  assembled  by  associating  a  carbon 
composite  electrode  with  a  MnC>2  electrode  in  order  to  increase 
the  cell  voltage  to  values  as  high  as  2.2  V.  Due  to  its  potential 
window  of  electroactivity,  the  MnC>2  electrode  is  used  as  the 
positive  electrode  and  carbon  as  the  negative. 

For  clarity  purpose,  the  design  shown  in  Fig.  1  mentions  two 
activated  carbon  electrodes.  However,  models  and  simulations 
developed  in  this  study  can  be  applied  to  any  type  of  cell,  such 
as  hybrid  aqueous  activated  carbon/Mn02  supercapacitor  devel¬ 
oped  by  our  team  (Tabema  et  al.  [7]).  Thermal  properties  of 
negative  and  positive  electrode  depend  on  their  structures,  espe¬ 
cially  the  influence  of  polymer  and  carbon  contents. 

Dimensions  and  properties  of  the  various  components  are 
given  in  Table  1 . 

4.  Finite  element  analysis 

For  optimal  design  purposes,  it  is  desirable  to  know  how  the 
temperature  of  an  electrochemical  capacitor  varies  in  space  and 
time  in  any  point  of  interconnected  cells.  The  relation  between 
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Table  1 


Characteristics  of  carbon-carbon  and  carbon-Mn02  supercapacitors,  respec¬ 
tively,  from  [7,8] 


Activated  car¬ 

bon/activated 

carbon 

Activated 

carbon/Mn02 

Specific  capacitance 

95  Fg-1  (AC) 

150  Fg-1  (Mn02) 

(positive  electrode) 

Specific  capacitance 

95  Fg"1  (AC) 

105  Fg-1  (AC) 

(negative  electrode) 

Electrolyte 

N(et)4+BF4  in 
acetonitrile 

K2SO4  in  water 

Cell  voltage  (V) 

2.3 

2.2 

Equivalent  serie  resis- 

^1 

~2 

tance  (Q  cm2) 

Power  ( kW  kg~ 1 ) 

44 

30 

Energy  (Whkg-1) 

17 

15 

Charge-discharge  test 

<1 

<1 

current  (A  cm  2) 

the  heat  production  rate  during  repetitive  charge-discharge  and 
the  temperature  T(M,  t)  at  any  point  of  the  electrochemical  cells 
is  basically  ruled  by  the  heat  diffusion  equation: 


Prescribed  temperature :  25°C 

Max:  33.527 


1 - ► 

O  x  (m) 


Fig.  2.  Temperature  distribution  in  the  unit  cell,  per  1 W  per  cell  heat  source. 

4.1.  Temperature  distribution  of  the  electrochemical  unit 
cell 


V  •  (— A(m)  •  Vf(M,p)  +  cp — =  Qy 

ot 


The  heat  generation  rate  Q  is  caused  by: 


•  ionic  transport  in  electrolyte  (electrodes  and  separator)  and 
electronic  charge  transport  in  collectors  and  solid  phase  of 
electrodes; 

•  reversible-irreversible  electrochemical  reactions  at 
solid-liquid  interface  of  the  porous  structures; 

•  thermal  contact  and  electrical  resistances  between  layers. 


Electrochemical  reactions  occur  in  the  positive  electrode  of 
the  hybrid  cells. 

During  periodic  charge-discharge  current  cycles,  the  heat 
generation  rate  can  be  decomposed  into  two  components: 

G(0  =  Q  +  8Q(t)  (2) 

where  Q  is  the  mean  value  reached  in  steady- state  conditions  for 
the  periodic  repetition  rate  of  the  current  cycle. 

Consequently,  the  temperature  at  any  point  of  an  electro¬ 
chemical  capacitor  can  be  written  as: 

T  =  T  +  8T(t )  (3) 


An  interconnected  bipolar  electrochemical  cell  is  considered 
in  this  study,  in  the  case  where  the  current  flows  normal  to  the 
current  collector  from  one  cell  to  the  next.  Temperature  is  cal¬ 
culated  at  any  point  inside  the  cell  using  the  steady- state  heat 
conduction  equation: 

-A,  V2r  =  Qi  (4) 

Five  different  layers  are  considered:  half  of  the  two  metal  collec¬ 
tors,  the  negative  and  positive  electrode  and  the  separator.  The 
cell  is  assumed  to  be  thermally  insulated,  except  the  two  terminal 
ends  of  the  collectors  (Fig.  2).  The  temperature  of  the  two  termi¬ 
nal  ends,  the  so-called  case  temperature,  is  prescribed.  Electrical 
energy  is  stored  in  a  double  layer  at  the  solid-electrolyte  inter¬ 
face  of  the  porous  electrodes.  Heat  dissipation  is  produced  into 
the  volume  of  individual  cells,  porous  electrodes,  electrolyte  and 
collectors  (see  Table  2).  The  two-dimensional  temperature  distri¬ 
bution  was  obtained  by  the  finite  element  method  implemented 
in  Femlab®/Matlab®  environment. 

Because  of  the  extreme  slimness  of  the  cell,  the  temperature 
gradient  is  found  to  be  very  weak  along  the  Ox  axis  (Fig.  2).  The 
huge  thermal  transfer  surface  between  two  adjacent  layers  favors 
heat  conduction  in  the  Ox  direction  and  allows  the  Oy  gradient 
to  be  high.  The  maximum  temperature  is  located  at  the  medium 


where  T  is  the  steady-state  temperature.  In  most 
charge-discharge  repetition  rates,  dynamic  temperature 
variations  can  be  neglected  because  the  heat  diffusion  time 
constant  (r^  10  min  in  the  case  described  in  this  paper)  is 
largely  greater  than  the  electrical  pulse  period. 

The  paper  investigates  the  temperature  distribution  T  at  any 
point  of  electrochemical  cells  considering  Q  as  a  known  space 
function.  Q  is  related  to  charge  transport  phenomena  (see,  for 
example,  Guillemet  et  al.  [8]). 


Table  2 


Thermophysical  properties  of  the  electrochemical  unit  cell 


Thermal  conductivity 
MlWm-'r1) 

Heat  source  (W) 

Collectors 

200 

0.05  each 

Electrodes 

0.5 

0.20  each 

Separator 

0.9 

0.50 

Total 

1 

Total  dissipated  power  is  1 W. 
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Fig.  3.  Temperature  inside  the  cell,  along  collector  and  separator,  per  1  W  per 
cell  heat  source. 

line  of  the  cell  (y  =  lyl 2  =  4.5  cm).  The  hot  spot  temperature  is 
quite  important  (8.5  K  W_1  above  the  case  temperature)  in  the 
conditions  of  the  simulation  (Fig.  3).  Temperature  is  roughly  the 
same  along  each  layer  of  the  cell. 

4.2.  Association  of  cells 

In  order  to  attain  high  values  of  capacitance,  elementary  cells 
have  to  be  stacked,  typically  a  few  hundreds  cell  associations. 

In  order  to  investigate  the  thermal  behavior  of  such  structures, 
a  finite  element  method  (FEM)  computation  is  carried  out  within 
an  association  of  50  cells,  with  an  electrical  insulator  (thick¬ 
ness:  ^is  =400  (xm,  thermal  conductivity:  =  0.4  Wm_1  K_1) 
and  an  external  package  (thickness:  eb  =  2mm,  thermal  con¬ 
ductivity:  Ab  =  200  Wm-1  K-1).  The  temperature  is  prescribed 
at  the  electrical-thermal  interconnects  and  the  external  surface 
(T0  =  25  °C). 

The  FEM  computation  in  Fig.  4  shows  that  the  temperature 
distribution  of  the  central  cell  is  not  really  different  from  the  one 
investigated  in  the  previous  section.  However,  the  cell  located  at 
the  edge  is  greatly  influenced  by  the  thermal  conduction  through 
the  package.  The  high  thermal  conductivity  of  the  package  leads 
to  a  quasi-uniform  temperature  in  the  near  cell,  as  expected. 

Electrical-thermal  interconnects 


Fig.  4.  Temperature  distribution  in  the  association  of  50  unit  cells  (FEM  com¬ 
putation,  per  1  W  per  cell  heat  source). 


The  maximum  temperature  is  reached  at  the  center  of  the  cell 
association.  The  hot  spot  temperature  is  8.4  °C  over  the  case 
temperature,  which  is  very  close  to  the  over  heating  computed 
in  an  elementary  cell  in  the  previous  section. 

5.  Effect  of  external  cooling  conditions 

External  convection-radiation  cooling  is  commonly  consid¬ 
ered  as  the  way  of  avoiding  internal  overheating  of  the  cells. 
Such  a  cooling  mode  basically  consists  of  air  flow  at  ambient 
temperature.  The  cooling  rate  is  then  related  to  the  local  heat 
transfer  coefficients  at  the  surface  package  and  terminal  con¬ 
nections.  The  heat  transfer  coefficient  on  a  plane  wall  in  the 
air  at  25  °C,  denoted  by  /zb,  is  typically  between  5  W  m-2  K-1 
for  natural  convection  and  200  W  m~2  K_1  for  turbulent  forced 
convection. 

The  effect  of  thermal-electrical  interconnects  with  the  supply 
wires  is  often  underestimated.  In  the  case  of  an  ultra  capaci¬ 
tor,  the  large  contact  surface  between  electrodes  and  collectors 
and  the  high  thermal  conductivity  of  the  collector  makes  the 
collectors  a  very  efficient  heat  spreader  (Fig.  5)  only  if  the 
electrical-thermal  interconnects  are  especially  designed. 

5.1.  Terminal  conductance  of  electrochemical  capacitors 

Considering  the  equivalent  thermal  circuit  of  the  geometry  of 
Fig.  5,  a  thermal  resistance  Rj  between  end  collectors  and  sur¬ 
rounding  air,  can  be  calculated  neglecting  thermal  resistance  of 
metal  bounding  and  thermal  contact  resistance  between  bound¬ 
ing  and  collectors: 

Rj  =  Rc  +  Rw  (5) 

where  Rc  is  the  thermal  resistance  of  N  collector  cells  connected 
in  parallel  to  the  metal  bounding: 

6c  1 

R<  =  VTT  x  (6) 

hclcdc  14 

and  where  Rw  is  the  equivalent  thermal  resistance  of  the  terminal 
wire,  subjected  to  air  flow,  calculated  from  the  basic  approxima¬ 
tion  of  extended  surfaces  [9] .  The  terminal  wire  is  considered  as 


Fig.  5.  Thermal  drain  due  to  electrical  connections  (electrical  connections  on 
lower  side  of  cells  are  identical  and  not  drawn). 
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an  infinite  fin  of  uniform  cross-sectional  area: 


1 


KM 


The  fin  factor  M  is: 


M  = 


w^w 


Using  common  values:  hw  =  100  Wm-2  K_1,  K  = 
200  W  K1  m1 ,  dw  =  5  cm,  ec  =  60  (Jim,  ac  =  200  W  K1  m  1 , 
/C  =  3cm,  dc  =  1  cm,  N=  50,  one  obtains:  Rtot  =  0.96KW_1. 
The  equivalent  conductance  can  be  considered  as  a  collector 
boundary  conductance  (or  terminal  conductance  Kt)  in  the 
thermal  model,  uniformly  distributed  on  the  whole  upper  (or 
lower)  surface  of  the  collectors: 


R 


tot 


N  X  ( ly(?c ) 


The  order  of  magnitude  of  Ky  is  thus  3860  Wm-2  K-1.  This 
is  quite  a  high  conductance  value,  much  higher  than  ones  with 
convection  flow  on  a  simple  wall.  The  origin  of  these  phenomena 
is  due  to  the  large  exchange  surface  offered  by  the  terminal 
wire  in  air.  Therefore,  in  the  case  of  natural  convection  cooling, 
the  terminal  conductance  value  can  be  estimated  around  few 
100  W  m  2  K-1,  and  in  case  of  turbulent  convection  cooling,  it 
might  be  possible  to  reach  20,000  W  m-2  K_1 . 


5.2.  Hot  spot  temperature  as  a  function  of  external  cooling 
conditions 

In  the  FEM  computations  above,  a  prescribed  temperature 
has  been  applied  on  the  terminal  end  of  collectors  as  a  boundary 
condition.  It  means  an  infinite  value  of  Kj  which  is  not  reachable 
in  practical  situations.  Using  different  values  of  terminal  con¬ 
ductance  and  heat  transfer  coefficient  in  a  third  kind  of  boundary 
condition  on  the  collectors,  different  hot  spot  temperatures  are 
obtained,  in  Figs.  6  and  7. 

Fig.  6  shows  that  when  the  electrical  connections  are  ther¬ 
mally  improved  (highest  terminal  conductance  Kj ),  the  hot  spot 
temperature  is  lower,  as  expected,  and  mainly  independent  of 
the  external  cooling  rate  (h^  value):  most  of  heat  is  evacuated 
through  electrical  terminal  ends. 

Fig.  7  shows  that  for  any  value  of  the  hot  spot  temper¬ 
ature  remains  strongly  dependent  on  the  value  of  Kt .  In  those 
conditions,  it  is  obvious  that  a  poor  electrical  contact  between 
the  terminal  wire  and  the  collector  ends  not  only  generates  an 
additional  heat  source,  but  also  adds  a  thermal  resistance  that 
lowest  value  of  Kt . 

This  asserts  that  great  attention  has  to  be  given  to  thermal 
design  of  the  electrical  connections  as  well  as  to  their  cooling. 


Fig.  6.  Maximum  temperature  in  the  cell  collection  vs.  boundary  cooling  of 
package,  for  extreme  boundary  conditions  of  collector,  per  1  W  per  cell  heat 
source. 


hB(W.m'2.K'') 


Fig.  7.  Maximum  temperature  in  the  cell  collection  vs.  boundary  cooling  of 
collectors,  for  extreme  boundary  conditions  of  package,  per  1  W  per  cell  heat 
source. 


6.  Correlation  between  surface  and  hot  spot 
temperature 

Considering  that  the  temperature  dependence  of  super  capac¬ 
itor  parameters  (capacitance,  resistance,  . . .)  requires  that  the 
definition  of  the  “cell  temperature”  is  clear  leads  to  the  problem 
of  the  temperature  measurement.  The  most  easy  and  most  used 
way  to  measure  the  “super  capacitor  temperature”  is  to  tape  a 
thermocouple  or  a  platinum  probe  to  an  electrical  connection  of 
the  package,  close  to  the  collectors.  Because  the  thermal  resis¬ 
tance  between  the  center  of  each  cell  and  the  end  of  the  collector 
is  not  negligible,  the  maximum  cell  temperature  is  not  the  one 
that  is  measured  (see  Fig.  8).  Moreover,  the  difference  between 
those  two  temperatures  cannot  be  considered  as  constant  in  so  far 
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Fig.  8.  Error  measurement,  between  hot  spot  temperature  and  collector  temper¬ 
ature  vs.  collector  heat  transfer  coefficient,  per  1  W  per  cell  heat  source. 

as  it  depends  widely  on  the  value  of  the  heat  transfer  coefficient 
/zb  and  conductance  Kj. 

Fig.  8  points  out  that  the  error  measurement  (Tmax  —  7c) 
is  weak  and  independent  of  the  cooling  rate  of  terminal  end 
of  the  collectors  ( Kj )  when  the  cooling  rate  of  the  package  is 
weak  (/zb  =  5  W  m-1  K_1),  but  increases  roughly  and  becomes 
strongly  dependent  of  Kj  in  the  case  of  an  efficient  package 
cooling  rate  (/zb  =  100  Wm-1  K_1). 

7.  Reduced-order  thermal  modeling 

Besides  FEM  computations,  there  is  a  need  for  a  simple  but 
accurate  thermal  model  that  could  be  used  in  the  design  of  ultra 
capacitors  and  for  prediction  of  hot  spots. 

7.7.  Heat  sources  equivalent  distribution 

Because  of  the  weak  Ox  temperature  gradient,  the  tempera¬ 
ture  distribution  of  the  cell  can  be  calculated  under  the  assump¬ 
tion  that  heat  sources  are  localized  inside  the  collectors.  To 
validate  this  assumption,  a  simulation  was  done,  close  to  the 
simulation  described  in  Section  3. 

The  total  heat  source  is  kept  (1  W  per  cell),  but  this  power  den¬ 
sity  is  no  longer  distributed  in  the  separator  and  the  electrodes, 
but  concentrated  in  the  collectors  (see  Table  3).  Results  show  that 
this  approximation  leads  to  an  overestimation  of  only  0.2%  of  the 
hot  spot  temperature  rise.  The  good  agreement  between  the  exact 
distribution  of  the  heat  source  and  localized  heat  source  approx¬ 
imation  reduces  the  problem  of  identification  of  microscopic 

Table  3 

Power  density  repartition  in  the  layers  of  the  cell,  in  case  of  distributed  heat 
sources  or  localized  heat  source  approximation 


Distributed  sources 

Localized  sources 

Collectors  (W  m-3) 

2.00  x  105 

2.00  x  106 

Electrodes  (W  m-3) 

1.85  x  105 

0 

Separator  (W  m~3) 

1.20  x  106 

0 

Total  (W) 

1 

1 

T max  —  Timh  (  C) 

8.5 

8.5 

hB=100  W.m'2.K' 
hB=5  W.m'lK/1 


1000 


10000 

K.T(W.m'2.K.') 


100000 


Fig.  9.  Description  of  the  association  of  elementary  cells  as  a  shell  network. 

inner  heat  sources  due  to  volume  joule  effects,  thermal  contact 
resistance,  electrochemical  reactions, ....  A  global  calorimetric 
measurement  is  sufficient,  a  priori,  to  characterize  all  of  those 
heat  sources. 

7.2.  Shell  network  modeling 

In  this  model,  the  association  of  A  elementary  cells  is  consid¬ 
ered,  in  steady-state  conditions,  as  an  array  of  Acollectors  where 
total  heat  generation  rate  is  concentrated,  separated  by  an  equiv¬ 
alent  resistive  element  that  models  the  two  electrodes  and  the 
separator  (Fig.  9).  Therefore,  the  equivalent  thermal  conductiv¬ 
ity  of  the  multi-layered  structure  electrode-separator-electrode 
in  the  Ox  direction  from  one  collector  to  the  next  is  introduced. 
The  equivalent  thermal  circuit  leads  to  the  following  relation: 


^es1  —  2/3  •  7e  1  +  (1  —  2/3)  •  7S  1 


(10) 


with  /3  = 


f  e+f  s 


All  of  the  elements  are  considered  as  a  homoge¬ 


neous  medium.  The  boundary  conditions  are:  prescribed  temper¬ 
ature  on  one  of  the  two  extremities  of  each  collector  (or  infinite 
value  of  hc)  and  on  metal  package  side  walls,  thermal  insulation 
elsewhere.  These  boundary  conditions  are  the  best  to  test  the 
validity  of  the  thermal  models  because  they  involve  the  highest 
temperature  gradients  in  the  medium. 

An  energy  balance,  based  on  the  extended  surface  approxi¬ 
mation,  so-called  shell  approximation  [9],  leads  to,  for  the  first 
element  (surface  package): 


d27 Ky) 
dy2 


-  (Kb  +  yo)Ti(y)  +  K0?2(y)  =  -n^mb 


(11) 


with  y\,  =  fj-  and  yo  =  yrfj-  •  Then,  for  the  following  elements 
( N  collectors),  from  z  =  1  to  N : 


dzTi(y )  Q 

- -f1  +  yTi-i(y)  -  2yTi(y)  +  yTi+x(y)  = 


dyz  '  '  7C 

with  y  =  and  for  the  last  element  (surface  package): 

d27yv+i(y) 


(12) 


dy: 


+  VoTN(y)  -  (Kb  +  Ko)Ta'+i(k)  =  -KbTamb  (13) 
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The  simultaneous  equations  (Eqs.  (1 1)— (13))  can  be  written  as: 
d2|T(y)] 

,  /  -  [G][roo]  =  [Q]  (14) 

dyz 

where  [7]  is  the  temperature  vector  Tfy).  The  solving  method 
uses  eigenvalue  matrix  computation.  In  the  eigenvectors  space, 
Eq.  (14)  becomes: 

d2[6>0(y)] 

.  ,  -  [Go][0o(j)]  =  [Qo]  (15) 

dyz 

where  [Go]  is  the  eigenvalues  diagonal  matrix,  rjj ,  of  the  matrix 
[G]  and  [R]  is  the  transformation  matrix,  as: 

[G0]  =  [/5]“1[G][P]  (16) 

The  matrix  [$ofe)]  and  [Qo]  are: 

[Qo]  =  [P]~1lQ]  and  [0o(z)]  =  m_1  [roo] 

The  general  solution  of  Eq.  (15)  is,  from  i  =  0  to  (N+  1): 

Ooiiy)  =  Oich((y  -  L)mi )  +  Aish((y  -  L)nti )  -  (17) 

mf 

where 


Table  4 

Recap  of  the  different  “hot  spot  to  case”  thermal  resistances  obtained  by  the 
different  models 


Thermal  resistance  (K  W  1 ) 

Finite  element  method 

8.4 

Localized  sources 

8.4 

Shell  network  model 

8.2 

Homogenization  method 

8.2 

Ultra-reduced-order  model 

7.2 

The  thermal  power  is  considered  to  be  1  W  per  cell  and  can  reach  few  Watts  per 
cell  in  practical  situations. 


The  thermal  conductivity  in  the  transverse  direction  (y  axis) 
is  then: 

Xy  =  ot  •  Xc  +  (1  —  ot)  •  Xes  (22) 

Using  the  values  given  in  Section  3  (see  Table  4),  one  obtains: 
AJC  =  0.64Wm_1  K"1  and  Ay  =  32.9  W  m- 1  KG1. 

An  FEM  computation  of  the  heat  conduction  equation  in 
the  cell  with  equivalent  Ox  and  Oy  direction  thermal  conduc¬ 
tivities,  thermal  insulation  on  side  boundaries  and  prescribed 
temperature  on  upper  and  lower  boundaries  (Fig.  10)  leads  to 
the  temperature  field  shown  in  Fig.  11. 


mi  = 


(18) 


The  constant  factors  Q{  and  A[  are  determined  from  boundary 
conditions: 


T  (0)  =  Tamb  and 


my) 

d  y 


=  0 


y=L 


for  even  values  of  i,  and: 


T(L)  —  Tamb  and 


d  T(y) 

d y 


=  0 


y=0 


(19) 


(20) 


for  odd  values  of  i. 

The  temperature  can  be  calculated,  then,  in  each  collector, 
along  the  Oy  axis.  Results  show  a  good  agreement  with  the 
FEM  computation,  except  a  moderate  underestimation  of  the  hot 
spot  temperature  (7max  —  Tamb  =  8.2  °C  for  shell  network  model, 
8.4  °C  for  FEM  computation). 


7.3.  Homogenization  method 

It  is  of  interest  to  apply  the  standard  homogenization  method 
in  the  case  of  the  ultra  capacitor  layout.  This  method  is  com¬ 
monly  used  for  the  thermal  modeling  of  various  electrochemical 
storage  systems  (see,  for  example,  refs.  [10,11]). 

The  multi-layered  structure  of  the  capacitor  is  replaced 
by  homogeneous  orthotropic  structure  having  equivalent  ther¬ 
mal  properties  and  where  the  heat  source  is  uniformly  spread 
throughout  the  whole  volume.  The  thermal  conductivity  in  the 
direction  normal  to  the  collectors  (x  axis)  is  then: 

Xx  =  ot  •  kc  T  (1  —  ot)  •  A.^1  (21) 

with  a  =  and  where  Aes  is  defined  as  in  the  previous  sec¬ 

tion. 
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Fig.  10.  Conductivities  and  geometry  of  an  elementary  cell  (top)  and  equivalent 
homogeneous  element  (bottom). 


Ph.  Guillemet  et  al.  /  Journal  of  Power  Sources  157  (2006)  630-640 


637 


Fig.  1 1 .  Temperature  field  in  the  elementary  cell,  atx  =  lxl 2,  comparison  between 
the  exact  heterogeneous  model  and  the  homogeneous  approximation. 


Using  the  values  given  in  Section  4.2,  the  hot  spot  temperature  is 
now  8.2  °C  over  ambient  temperature,  instead  of  8.4  °C  obtained 
by  FEM  computation  on  a  heterogeneous  model,  applied  to  the 
association  of  Section  4.2.  This  result  is  quite  good:  the  homoge¬ 
nization  method  is  fully  applicable  to  the  case  of  ultra  capacitors. 

7.4.  Ultra-reduced-order  modeling 


The  maximum  temperature  elevation  is  found  to  be  10.0  °C 
in  the  homogeneous  approximation,  instead  of  10.6  °C  in  the 
exact  heterogeneous  model,  applied  to  a  unique  cell  (Section  3). 

In  the  case  of  the  association  described  in  Fig.  12,  the  cal¬ 
culation  is  nearly  the  same.  The  thermal  conductivities  of  the 
whole  homogeneous  medium  are  the  same  as  for  a  single  cell 
(kx  =  0.643  W  m~ 1  K- 1  and  ky  =  32.9  W  m~ 1  K~ 1 )  but  the  insu¬ 
lator  and  package  layers  have  to  be  taken  into  account  through 
thermal  resistances.  The  surface  conductance  is  then: 


1 

Ky  -  TjTTVT 

A-b  A-is 


(23) 


7.4.1.  Four  port  matrix  formalism  for  ID  heat  conduction 
analysis 

Four  port  matrix  formalism  can  be  used  for  ID  heat  conduc¬ 
tion  analysis  with  internal  heat  production. 

Consider  the  ID  heat  conduction  in  a  slab  of  homogeneous 
medium  of  Fig.  13.  Solution  of  the  heat  conduction  equation 
leads  to  Eqs.  (A. 8)  and  (A. 11),  giving  relations  between  tem¬ 
perature  and  heat  flux  at  terminal  ends  and  heat  generation  rate: 

^out  —  —  Qtot^  ^^in  T  ^in 

^out  =  Qtot  +  ^in 


Eqs.  (A. 8)  and  (A.  11)  can  be  transformed  in  a  matrix  relation, 
between  in  and  out  values: 


This  relation  describes  the  thermal  behavior  of  the  slab. 


7.4.2.  Equivalent  thermal  circuit 

Consider  the  thermal  circuit  of  Fig.  14.  The  energy  balance 
is  obviously: 

^out  —  Gtot  T  ^in  (25) 


which  is  no  more  than  Eq.  (A.l  1).  The  incoming  flux  is: 


Fig.  12.  A  few  cells  from  the  association  of  50  elementary  cells. 


R 

2 


(26) 
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O  T  O 

in  0  out 


dr 
—A  — 
d y 


y= o 


0 


Jin 


Ix^z 


KpiTyO  7air) 


dr 
—A  — 
d  y 


0 


Jout 


Ix^z 


KT(Ty0  rair) 


(30) 

(31) 


Fig.  14.  Equivalent  thermal  circuit  from  the  homogeneous  slab  of  Fig.  13. 


A  heat  flux  balance,  at  the  node  O ,  leads  to: 


Tin  To  Tou{  To 

j  +  ytot  +  j 

2  2 


(27) 


Eliminating  To  between  those  two  relations  leads  exactly  to  Eq. 
(A. 8):  the  circuit  of  Fig.  14  is  the  equivalent  thermal  circuit  of 
the  slab  in  Fig.  13,  i.e.  to  the  matrix  equation  (Eq.  (24)),  in  term 
of  input-output  quantities.  It  enables  the  calculation  of  in-  and 
out- temperature  and  fluxes  ( Tm ,  Tout,  0m ,  T>oui)  but  not  of  the 
temperature  inside  the  medium:  the  temperature  7o  at  the  central 
node  of  Fig.  14  is  not  a  physical  temperature  but  a  mathematical 
one,  on  a  virtual  point.  Once  the  input  and  output  quantities  have 
been  determined,  the  hot  spot  temperature  has  to  be  calculated 
using  the  heat  conduction  equation  in  the  medium,  which  can 
be  very  easy  as  shown  below. 


In  this  diagram,  R(  is  the  thermal  resistance  in  the  Ox  or  Oy 
direction: 


Rx  = 


Nl 


X 


TXlylZ 


(32) 


Ry  = 


l 


y 


TylXlZ 


(33) 


The  thermal  conductivities  Xx  and  Xy  have  been  calculated  in 
Eqs.  (21)  and  (22). 

Eq.  (36)  is  the  same  as  Eq.  (23),  it  takes  into  account  the  ther¬ 
mal  resistance  of  the  insulator  and  package  walls.  Calculation 
of  To  is  now  obvious,  using  the  energy  balance  at  the  node  O , 
for  example: 


Qtot  2 


To  T% mb 

Px  I  1 

2  hftlylz 


-2 


To  T% mb 

Py  I  1 

2  KtUz 


=  0 


(34) 


7.4.3.  Application  to  an  ultra  capacitor 

The  reduced-order  model,  detailed  in  the  above  section,  is 
convenient  and  easy  to  use.  It  enables  one  to  get  quickly  a  rough 
approximation  of  the  hot  spot  temperature  in  an  ultra  capacitor. 

The  equivalent  thermal  circuit  of  the  association  described  in 
Fig.  4  is  shown  in  Fig.  15.  It  takes  into  account  the  heat  flow  in 
both  Ox  and  Oy  directions.  Four  additional  thermal  resistances 
have  been  included  in  the  equivalent  circuit.  They  represent  the 
boundary  conditions: 

=  TT=  hB(Tx o  -  rair)  (28) 

x=o  hT 


According  to  the  previous  boundary  conditions  (prescribed  tem¬ 
perature  on  upper  and  lower  sides  and  on  external  package 
walls),  the  two  surface  conductances  become: 


K,  =  oc  (35) 

hB  =  (36) 

A-b  A-is 

Solving  this  equation  leads  to  the  value  of  To  which  has  no 
physical  meaning  but  it  enables  one  to  calculate  Txo  and  Ty o 
(see  Fig.  15): 


x — lx 


0 


•^out 


lylZ 


h  B  (  Txo  Tair) 


(29) 


T 


Aim  b  Txo  — 


Tamb  TyO  — 


Tamb  To 


2  hB^ylz 
Tamb  To 


Ry  1 

2  Kplxh 


(37) 


Using  again  the  values  given  in  Section  4.2,  the  temperatures 
are:  Txo  =  25.9  °C  and  Ty o  =  25.0  °C  (prescribed  temperature). 

To  complete  the  calculation,  the  hot  spot  temperature,  in  the 
middle  of  the  cell  association,  can  be  evaluated  with  a  basic 
approximation.  The  application  of  Eq.  (A. 4),  in  the  ID  heat 

transfer  it  describes,  with  v  =  lxl 2  leads  to  Qi0 1— 4  X=l/1R — —  =  0. 

2 

In  the  2D  heat  transfer  of  Fig.  15,  this  expression  can  be  extended 
as: 


Qtot  4 


Tmax 


Rx 

2 


Tmax  T 


yO 


R y 


(38) 


Fig.  15.  Equivalent  thermal  circuit  from  the  ultra  capacitor  of  Fig.  4. 
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where  Tmax  is  the  hot  spot  temperature.  Using  again  the  val¬ 
ues  given  in  Section  4.2,  one  obtains:  rmax  —  Tam b  =  7.2  °C.  The 
discrepancy  between  this  result  and  the  FEM  result  (8.4  °C)  is 
about  16%  which  is  not  too  much,  considering  the  ease  of  this 
method  to  calculate  the  hot  spot  temperature  using  only  three 
basic  equations  (Eqs.  (34),  (37)  and  (38)). 

8.  Summary  of  models  and  results 

The  numerical  FEM  computation  of  Section  4.2  and  the 
four  semi-analytical  models  developed  in  Section  7  have  both 
been  conducted  in  the  same  configuration  of  A  =50  associ¬ 
ated  cells  (Fig.  12),  with  a  heat  source  of  1  W  per  cell  and 
with  a  prescribed  connection  temperature,  involving  the  most 
important  temperature  gradients  in  the  Oy  and  Ox  directions 
and  leading  to  the  determination  of  the  “hot  spot  to  case” 
resistance. 

The  finite  element  method  gives  the  most  accurate  results  of 
the  temperature  distribution  anywhere  in  the  medium,  but  it  is  a 
weighty  numerical  method. 

The  shell  network  model  enables  the  analytical  computation 
of  the  temperature  distribution  in  the  collectors  and  leads  to 
a  15%  underestimation  of  the  hot  spot  temperature  rise  with 
respect  to  ambient  temperature. 

The  homogenization  method  gives  very  good  results  (less 
than  2%  underestimation  in  the  studied  case).  Although  it  is  a 
numerical  method,  this  approximation  leads  to  a  large  decrease 
of  computing  time  compared  to  the  FEM  method,  considered  as 
reference  in  a  lack  of  experimental  measurements. 

The  ultra-reduced-order  model  offers  the  great  advantage  of 
providing  an  easy  and  fast  analytical  calculation  of  the  hot  spot 
temperature  in  an  ultra  capacitor:  this  model  is  to  be  used  in 
most  practical  case  where  a  high  accuracy  is  not  needed.  For 
example,  thermal  security  predictions  always  include  a  large 
safety  range  which  is  usually  wider  than  the  accuracy  of  the 
model. 

9.  Application  to  a  parametric  study 

The  reduced- order  model  leads  to  a  rapid  and  accurate  esti¬ 
mation  of  the  hot  spot  cell  temperature,  depending  on  a  wide 
range  of  parameters:  the  cooling  rate,  the  properties  and  size  of 
constituting  materials,  the  dissipated  power  source  as  function 
of  repetition  rate  of  current  cycles. 

It  is  of  interest,  for  example,  to  evaluate  the  influence  of  col¬ 
lector  conductivity  on  this  temperature.  Eqs.  (21)  and  (22)  give 
the  dependence  of  the  conductivity  values  of  the  homogeneous 
equivalent  medium  Xx  and  Xy ,  then  Eqs.  (32)  and  (33)  give  the 
resistance  values  Rx  and  Ry  and  finally  Eq.  (38)  gives  the  hot 
spot  overheating  for  a  1  W  heat  source,  which  is  the  thermal 
resistance  of  the  ultra  capacitor. 

Fig.  16  shows  that  decreasing  the  collector  conductivity  from 
400  Wm-1  K-1  down  to  ^Wm"1  K_1  increases  the  thermal 
resistance  from  4  K  W- 1  up  to  22  K  W- 1 . 

Moreover,  the  comparison  between  this  basic  model  (so- 
called  ultra-reduced-order  model)  based  on  multi-port  element 
network  and  the  exact  FEM  computation  (Fig.  16)  shows  the 


Fig.  16.  Influence  of  collector  thermal  conductivity  on  the  “hot  spot  to  case” 
thermal  resistance. 

validity  of  the  approximation  over  a  wide  range  of  materials  and 
sizes. 

10.  Application  to  thermal  safety 

The  thermal  resistances  given  in  the  previous  section  were 
obtained  with  Q  =  1  W  per  cell  (8 1  cm2).  In  practice,  the  thermal 
power  can  be  much  more  than  1  W,  especially  in  high  inten¬ 
sity  current  applications  that  are  the  main  application  field  of 
ultra  capacitors.  Consider  the  case  of  an  ultra  capacitor  with 
5  £2  cm-2  ESR  and  0.33  F  cm-2  capacitance.  A  200  mA  current 
intensity  (equivalent  to  166  A  in  a  2500  F  association)  leads  to 
1 .8  W  per  unit  cell  and  consequently  to  about  15  K  overheating. 
This  overheating  can  be  detrimental  in  high  ambient  temperature 
(Tamb  =  40  °C)  and  low  cooling  rate  (rcase  =  60  °C).  In  those  con¬ 
ditions,  the  hot  spot  temperature  reaches  75  °C  and  can  overtake 
the  safety  limit  of  running  temperature  (70  °C  for  acetonitrile 
solvent)  and  cause  serious  damage. 

11.  Future  works 

The  present  work  is  a  preliminary  study  to  future  works  that 
intend  to  achieve  a  complete  description  of  heat  transport  phe¬ 
nomena  in  ultra  capacitors,  including  other  geometries,  other 
materials,  thermal  transient  analysis  as  well  as  thermal  measure¬ 
ment  (temperature  change,  dissipated  power,  thermophysical 
properties  of  constituting  layers). 

12.  Conclusion 

Thermal  analysis  of  super  capacitors  points  out  that  the  maxi¬ 
mum  temperature  is  reached  at  the  center  of  the  cell  association, 
as  expected.  This  maximum  temperature  is,  to  a  large  extent, 
different  from  the  collector  temperature  that  can  be  measured 
outside  the  package.  Moreover,  this  hot  spot  temperature  cannot 
be  easily  related  to  the  outside  temperature  because  it  depends 
widely  on  the  way  the  heat  is  thrown  away  from  the  cell  associ¬ 
ation. 

In  order  to  evacuate  heat  produced  from  cells,  special  atten¬ 
tion  must  be  given  to  the  cooling  performance  of  the  electrical 
connections  of  the  collectors,  where  the  most  important  part  of 
temperature  rise  is  located. 
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Temperature  inside  the  ultra  capacitor,  or  thermal  resistance, 
can  be  calculated  using  four  models,  corresponding  to  four  levels 
of  thermal  analysis: 

•  Finite  element  method  gives  the  most  accurate  temperature 
and  flux  values,  but  is  a  weighty  method. 

•  The  shell  network  model  leads  to  a  semi-analytical  calculation 
of  the  temperature  in  each  layer  of  the  capacitor  under  the 
assumption  that  heat  transfer  is  longitudinal  from  one  to  the 
next  collector  and  in  the  axial  direction  through  the  collectors. 

•  The  homogenization  technique  gives  a  good  estimation  of  the 
temperature  everywhere  in  the  capacitor,  assuming  that  the 
cell  association  is  a  homogeneous  non-isotropic  structure. 

•  The  ultra-reduced-order  model,  based  on  the  thermal  proper¬ 
ties  homogenization  and  an  equivalent  thermal  circuit  leads  to 
a  simple  but  rather  good  estimation  of  the  thermal  resistance 
of  the  capacitor  and  enables  easy  study  of  any  parameter  influ¬ 
ence  on  this  thermal  resistance.  This  model  is  to  be  used  in 
most  practical  situation  where  a  high  accuracy  is  not  needed, 
like  thermal  security  prediction,  for  example. 


The  integration  constants,  A  and  B ,  can  be  expressed  from  Eq. 
(A. 2): 
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Applying  Eq.  (A. 6)  with  x  —  lx,  one  obtains: 

Tout  =  ~Qtot~  B0[n  +  7in 

where  R  is  a  thermal  resistance: 
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and  Qtot  is  the  total  heat  power  (W)  in  the  slab: 

Qtot  =  Qlxlylz 

Applying  Eq.  (A. 7)  with  x  =  lx,  one  obtains: 


(A.  10) 


Appendix  A 


The  ID  heat  conduction  equation,  in  the  slab,  is,  per  volume 
unit: 


d  2T 

dx2 


(A.l) 


where  Q  is  the  volumetric  heat  source  (W  m  3).  In  a  general 
case,  among  the  four  boundary  equations: 


T(x  =  0)  =  7-n, 
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two  equations  constitute  the  boundary  conditions,  used  to  solve 
this  ID  problem. 

The  general  form  of  the  solution  of  Eq.  (A.l)  is: 


Q_ 

2X 
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^out  —  Gtot  H"  ^in  (A.  11) 

which  is  no  more  than  the  energy  balance  into  the  slab. 
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